# 敏感性分析.py - 优化版本（仅保留3PG模型计算功能）
from flask import Flask, Blueprint, render_template, request, jsonify
import json
import numpy as np
from SALib.analyze import morris
from SALib.sample.morris import sample
from skopt import gp_minimize
from skopt.space import Real
import time
from three_pg_model import ThreePGModel, SiteData, default_parameters
import logging
from datetime import datetime
from models import ProjectsModel, SiteModel, SpeciesModel
from exts import db

calc = Blueprint('sens', __name__, url_prefix="/sens")

# 全局变量和配置
CONFIG_PATH = "Sensitivity/参数范围.json"
OUTPUT_VARS_PATH = "Sensitivity/输出可选变量.txt"

# 默认输出变量列表
DEFAULT_OUTPUT_VARS = [
    "StemNo", "WF", "WR", "WS", "StandVol", "LAI", "avDBH", "BasArea",
    "Height", "MAI", "TotalW", "NPP", "GPP", "Transp", "ET", "ASW"
]

# 结果描述字典 - 用于添加中文描述
RESULT_DESCRIPTIONS = {
    "StemNo": "林木数量 (株/ha)",
    "WF": "叶生物量 (t/ha)",
    "WR": "根生物量 (t/ha)",
    "WS": "茎生物量 (t/ha)",
    "WL": "枯落物 (t/ha)",
    "StandVol": "林分材积 (m³/ha)",
    "LAI": "叶面积指数 (m²/m²)",
    "avDBH": "平均胸径 (cm)",
    "Irrig": "灌溉量 (mm)",
    "incrWS": "茎生物量增量 (t/ha)",
    "mortality": "死亡率 (%)",
    "selfThin": "自疏 (株/ha)",
    "meantemp": "平均温度 (°C)",
    "FR": "肥力等级",
    "CO2": "CO2浓度 (ppm)",
    "minASW": "最小可用土壤水 (mm)",
    "StandAge": "林龄 (年)",
    "BasArea": "断面积 (m²/ha)",
    "Height": "树高 (m)",
    "MAI": "平均年生长量 (m³/ha/年)",
    "MAIx": "最大平均年生长量 (m³/ha/年)",
    "ageMAIx": "最大MAI年龄 (年)",
    "SLA": "比叶面积 (m²/kg)",
    "CanCover": "郁闭度 (%)",
    "LAIx": "最大LAI (m²/m²)",
    "ageLAIx": "最大LAI年龄 (年)",
    "TotalW": "总生物量 (t/ha)",
    "AvStemMass": "平均单株茎重 (kg)",
    "Density": "密度 (株/ha)",
    "fracBB": "树皮比例",
    "fAge": "年龄修正因子",
    "fVPD": "VPD修正因子",
    "fTemp": "温度修正因子",
    "fCalpha": "CO2对量子效率影响",
    "fCg": "CO2对气孔导度影响",
    "fFrost": "霜冻修正因子",
    "fSW": "水分修正因子",
    "fNutr": "养分修正因子",
    "PhysMod": "生理修正因子",
    "GPP": "总初级生产力 (tC/ha)",
    "NPP": "净初级生产力 (tC/ha)",
    "RADint": "截获辐射 (MJ/m²)",
    "alphaC": "冠层量子效率",
    "Epsilon": "光能利用效率 (gDM/MJ)",
    "CVI": "冠层体积指数",
    "m": "生物量分配系数",
    "pR": "分配至根的比例",
    "pS": "分配至茎的比例",
    "pF": "分配至叶的比例",
    "pFS": "叶茎分配比",
    "gammaF": "叶凋落率 (1/月)",
    "wSmax": "最大单株茎重 (kg)",
    "gammaN": "死亡率 (%/年)",
    "supIrrig": "补充灌溉 (mm)",
    "RunOff": "径流 (mm)",
    "fRainInt": "雨水截留比例",
    "RainInt": "雨水截留量 (mm)",
    "CanCond": "冠层导度 (m/s)",
    "WUE": "水分利用效率 (gDM/kg水)",
    "ET": "蒸散发 (mm)",
    "Transp": "蒸腾 (mm)",
    "ASW": "可用土壤水 (mm)",
    "DayLength": "日长 (h)",
    "FrostDays": "霜冻日数 (天)",
    "SolarRad": "太阳辐射 (MJ/m²)",
    "Tav": "平均温度 (°C)",
    "VPD": "饱和水汽压差 (mBar)",
    "Rain": "降雨量 (mm)",
    "ltStemGR": "长期茎生长率 (t/ha/月)",
    "Litter": "凋落物 (t/ha)",
    "TotalLitter": "总凋落物 (t/ha)",
    "EpsilonStem": "茎生长效率 (tDM/ha/MJ)",
    "pooledSW": "汇集土壤水 (mm)",
    "fT": "温度函数"
}


def parse_bound(bound, params):
    """解析参数边界值，支持数值和引用"""
    if isinstance(bound, (int, float)):
        return float(bound)
    elif isinstance(bound, str):
        # 支持参数间引用（如MaxASW）
        ref_param = next((p for p in params if p["name"] == bound), None)
        if ref_param:
            return parse_bound(ref_param["upper_bound"], params)
        raise ValueError(f"参数 {bound} 不存在")
    else:
        raise ValueError("无效的边界值类型")


def load_parameters_config():
    """加载参数配置"""
    try:
        with open(CONFIG_PATH, "r", encoding="utf-8") as f:
            params_data = json.load(f)

        params = params_data["parameters"]
        # 两阶段解析：先处理数值型边界，再处理引用型边界
        for param in params:
            param["lower_bound"] = parse_bound(param["lower_bound"], params)
            param["upper_bound"] = parse_bound(param["upper_bound"], params)

        return params  # 只返回参数列表
    except Exception as e:
        print(f"加载配置文件失败: {str(e)}")
        return []


available_params = load_parameters_config()


def get_output_vars():
    """获取输出变量列表"""
    try:
        output_vars = []
        with open(OUTPUT_VARS_PATH, "r", encoding="utf-8") as f:
            for line in f:
                if 'value="' in line and not 'value=""' in line:
                    var = line.split('value="')[1].split('"')[0]
                    output_vars.append(var)
        return output_vars
    except Exception as e:
        print(f"读取输出变量列表失败: {str(e)}")
        return DEFAULT_OUTPUT_VARS


def create_default_site():
    """创建默认站点数据"""
    site = SiteData()
    site.Sites = "默认站点"
    site.Latitude = 0.0
    site.FR = 0.0
    site.Soil_type = "CL"
    site.Max_ASW = 200.0
    site.Min_ASW = 0.0
    site.Initial_stocking = 1000.0
    site.gammaNx = 0.0
    site.alphaCx = 0.06
    return site


def get_default_climate_data():
    """获取默认气候数据"""
    return {
        "TMAX1": [20.9, 21.21, 18.91, 15.79, 12.38, 10.24, 9.88, 11.12, 12.76, 15.32, 16.57, 18.5],
        "TMIN1": [8.86, 9.98, 8.39, 6.91, 5, 3.11, 2.91, 3.01, 4.23, 5.05, 6.7, 8.03],
        "PRECIP1": [71.05, 62.51, 73.28, 108.89, 125.83, 118.06, 149.54, 126.51, 121.72, 126.8, 107.8, 102.83],
        "EVAP1": [3.95, 3.9, 2.7, 1.61, 1.08, 0.61, 0.7, 1.07, 1.64, 2.32, 2.98, 3.56],
        "RADN1": [23.47, 19.74, 13.68, 8.48, 5.43, 3.87, 4.65, 7.51, 10.85, 16.07, 19.45, 21.16],
        "RAIND1": [13.32, 11.35, 14.51, 17.28, 19.04, 17.18, 20.31, 20.12, 19.63, 18.69, 17.59, 16.63],
        "FROST1": [0 for _ in range(12)]
    }

def get_special_climate_data(sitename):
    """根据sitename获取气候数据"""
    site = SiteModel.query.filter_by(sitename=sitename).first()
    if site:
        return {
            "TMAX1": [float(site.tmax1), float(site.tmax2), float(site.tmax3), float(site.tmax4), float(site.tmax5),
                      float(site.tmax6),
                      float(site.tmax7), float(site.tmax8), float(site.tmax9), float(site.tmax10), float(site.tmax11),
                      float(site.tmax12)],
            "TMIN1": [float(site.tmin1), float(site.tmin2), float(site.tmin3), float(site.tmin4), float(site.tmin5),
                      float(site.tmin6),
                      float(site.tmin7), float(site.tmin8), float(site.tmin9), float(site.tmin10), float(site.tmin11),
                      float(site.tmin12)],
            "PRECIP1": [float(site.precip1), float(site.precip2), float(site.precip3), float(site.precip4),
                        float(site.precip5), float(site.precip6),
                        float(site.precip7), float(site.precip8), float(site.precip9), float(site.precip10),
                        float(site.precip11), float(site.precip12)],
            "EVAP1": [float(site.evap1), float(site.evap2), float(site.evap3), float(site.evap4), float(site.evap5),
                      float(site.evap6),
                      float(site.evap7), float(site.evap8), float(site.evap9), float(site.evap10), float(site.evap11),
                      float(site.evap12)],
            "RADN1": [float(site.ra1), float(site.ra2), float(site.ra3), float(site.ra4), float(site.ra5),
                      float(site.ra6),
                      float(site.ra7), float(site.ra8), float(site.ra9), float(site.ra10), float(site.ra11),
                      float(site.ra12)],
            "RAIND1": [float(site.raind1), float(site.raind2), float(site.raind3), float(site.raind4),
                       float(site.raind5), float(site.raind6),
                       float(site.raind7), float(site.raind8), float(site.raind9), float(site.raind10),
                       float(site.raind11), float(site.raind12)],
            "FROST1": [float(site.frost1), float(site.frost2), float(site.frost3), float(site.frost4),
                       float(site.frost5), float(site.frost6),
                       float(site.frost7), float(site.frost8), float(site.frost9), float(site.frost10),
                       float(site.frost11), float(site.frost12)]
        }
    else:
        return get_default_climate_data()

def get_special_species_data(speciesname):
    """根据speciesname获取树种参数"""
    species = SpeciesModel.query.filter_by(speciesname=speciesname).first()
    if species:
        params = {
            "MaxAge": float(species.MaxAge),
            "nAge": float(species.nAge),
            "rAge": float(species.rAge),
            "fullCanAge": float(species.fullCanAge),
            "k": float(species.k),
            "gammaR": float(species.gammaR),
            "SWconst0": float(species.SWconst),
            "SWpower0": float(species.SWpower),
            "fCalpha700": float(species.fCalpha700),
            "fCg700": float(species.fCg700),
            "MaxIntcptn": float(species.MaxIntcptn),
            "LAImaxIntcptn": float(species.LAImaxIntcptn),
            "MinCond": float(species.MinCond),
            "MaxCond": float(species.MaxCond),
            "LAIgcx": float(species.LAIgcx),
            "BLcond": float(species.BLcond),
            "CoeffCond": float(species.CoeffCond),
            "y": float(species.Y),
            "Tmax": float(species.Tmax),
            "Tmin": float(species.Tmin),
            "Topt": float(species.Topt),
            "kF": float(species.kF),
            "pFS2": float(species.pFS2),
            "pFS20": float(species.pFS20),
            "aWs": float(species.aS),
            "nWs": float(species.nS),
            "pRx": float(species.pRx),
            "pRn": float(species.pRn),
            "m0": float(species.m0),
            "fN0": float(species.fN0),
            "fNn": float(species.fNn),
            "alphaCx": float(species.alpha),
            "poolFractn": 0.0,
            "gammaN1": float(species.gammaNx),
            "gammaN0": float(species.gammaN0),
            "tgammaN": float(species.tgammaN),
            "ngammaN": float(species.ngammaN),
            "wSx1000": float(species.wSx1000),
            "thinPower": float(species.thinPower),
            "mF": float(species.mF),
            "mR": float(species.mR),
            "mS": float(species.mS),
            "gammaF1": float(species.gammaFx),
            "gammaF0": float(species.gammaF0),
            "tgammaF": float(species.tgammaF),
            "SLA0": float(species.SLA0),
            "SLA1": float(species.SLA1),
            "tSLA": float(species.tSLA),
            "fracBB0": float(species.fracBB0),
            "fracBB1": float(species.fracBB1),
            "tBB": float(species.tBB),
            "rho0": float(species.rhoMin),
            "rho1": float(species.rhoMax),
            "tRho": float(species.tRho),
            "aH": float(species.aH),
            "nHB": float(species.nHB),
            "nHN": float(species.nHN),
            "aV": float(species.aV),
            "nVB": float(species.nVB),
            "nVN": float(species.nVN),
            "Qa": float(species.Qa),
            "Qb": float(species.Qb),
            "gDM_mol": float(species.gDM_mol),
            "molPAR_MJ": float(species.molPAR_MJ),
        }
        # 验证 fullCanAge
        if params["fullCanAge"] <= 0:
            print(f"警告：树种 {speciesname} 的 fullCanAge 无效 ({params['fullCanAge']})，使用默认值。")
            params["fullCanAge"] = default_parameters["fullCanAge"]  # 使用默认值 3
        # 验证 pFS2 和 pFS20
        if params["pFS20"] < params["pFS2"]:
            print(f"警告：树种 {speciesname} 的 pFS20 ({params['pFS20']}) 小于 pFS2 ({params['pFS2']})，使用默认值。")
            params["pFS2"] = default_parameters["pFS2"]  # 默认值 1
            params["pFS20"] = default_parameters["pFS20"]  # 默认值 0.15
        return params
    else:
        print(f"未找到树种 {speciesname}，使用默认参数。")
        return default_parameters

def create_3pg_model(params=None, site=None, climate_data=None, Benchmark_id=None):
    """创建3PG模型实例"""
    if Benchmark_id is None:
        if site is None:
            site = create_default_site()
        if climate_data is None:
            climate_data = get_default_climate_data()

        #logging.info(f"Creating model with params: {params}")
        return ThreePGModel(
            site_data=site,
            species="Trial",
            climate_data=climate_data,
            seedling_mass=1,
            month_planted="June",
            planted_year=1990,
            end_age=10,
            lookup=False,
            params=params if params else default_parameters,
            output_frequency="m"
        )

    else:
        Benchmark = ProjectsModel.query.get(Benchmark_id)
        modelInput = json.loads(Benchmark.data_input)

        print(modelInput)

        site = SiteData()
        site.Sites = modelInput["Site"]
        site.Latitude = modelInput["Latitude"]
        site.FR = modelInput["Fertility_rating"]
        site.Soil_type = modelInput["Soil_class"]
        site.Max_ASW = modelInput["Maximum_ASW"]
        site.Min_ASW = modelInput["Minimum_ASW"]
        site.Initial_stocking = modelInput["Initial_stocking"]

        species = modelInput["Species"]

        sitename = modelInput["Site"]
        climate_data = get_special_climate_data(sitename)

        seedling_mass = modelInput["Seedling_mass"]

        month_planted = modelInput["Month_planted"]
        planted_year = modelInput["Year_planted"]

        end_age = modelInput["End_age"]

        lookup = modelInput["Interpolate_lookups"]

        speciesname = modelInput["Species"]
        params = get_special_species_data(speciesname)

        return ThreePGModel(
            site_data=site,
            species=species,
            climate_data=climate_data,
            seedling_mass=seedling_mass,
            month_planted=month_planted,
            planted_year=planted_year,
            end_age=end_age,
            lookup=lookup,
            params=params,
            output_frequency="m"
        )

def process_model_results(model_results, output_vars):
    """处理3PG模型结果，仅保留最后一个时间点数据"""
    final_results = []
    for var in output_vars:
        if var in model_results and len(model_results[var]) > 0:
            value = model_results[var][-1]
            if isinstance(value, (np.int32, np.int64)):
                value = int(value)
            elif isinstance(value, (np.float32, np.float64)):
                value = float(value)
            readable_key = RESULT_DESCRIPTIONS.get(var, var)  # 中文描述
            final_results.append({
                "english": var,  # 英文变量名
                "chinese": readable_key,  # 中文描述
                "value": round(value, 4)  # 计算结果
            })
    return final_results


# 3PG模型敏感性分析路由
@calc.route('/analyze_3pg', methods=['POST'])
def analyze_3pg():
    try:
        data = request.get_json()
        parameters = data.get('parameters', [])

        # 提取参数名称和边界
        names = [p["name"] for p in parameters]
        bounds = []
        for p in parameters:
            param_config = next((item for item in available_params if item["name"] == p["name"]), None)
            if param_config:
                lb = float(param_config["lower_bound"])
                ub = float(param_config["upper_bound"])
                if lb >= ub:
                    error_msg = f"参数 {p['name']} 的边界非法: 下界 {lb} >= 上界 {ub}"
                    print(error_msg)
                    return jsonify({"error": error_msg}), 400
                bounds.append([lb, ub])
            else:
                bounds.append([-100, 100])

        # 定义问题
        problem = {
            'num_vars': len(parameters),
            'names': names,
            'bounds': bounds
        }
        print("生成3PG模型分析问题:", problem)

        # 获取输出变量
        output_vars = get_output_vars()
        print(f"分析的输出变量: {output_vars}")

        # 生成Morris样本
        X = sample(problem, N=100, num_levels=10)  # N=10 减少计算量，实际可调整

        # 对每个样本运行3PG模型
        Y_matrix = np.zeros((X.shape[0], len(output_vars)))  # 样本数 × 输出变量数
        for i, x in enumerate(X):
            # 将样本参数转换为3PG模型参数
            user_params = dict(default_parameters)
            for j, param_name in enumerate(names):
                user_params[param_name] = x[j]

            # 创建并运行3PG模型
            model = create_3pg_model(params=user_params)
            model_results = model.run3PG(output_vars=output_vars)

            # 提取每个输出变量的最后值
            for k, var in enumerate(output_vars):
                if var in model_results and len(model_results[var]) > 0:
                    Y_matrix[i, k] = model_results[var][-1]
                else:
                    Y_matrix[i, k] = 0  # 未返回结果时填充0

        # 对每个输出变量进行Morris分析
        mu_star_matrix = []
        sigma_matrix = []
        for k in range(len(output_vars)):
            Si = morris.analyze(problem, X, Y_matrix[:, k], print_to_console=False)
            mu_star_matrix.append(Si['mu_star'].tolist())
            sigma_matrix.append(Si['sigma'].tolist())

        print("3PG模型敏感性分析完成，矩阵大小:", len(mu_star_matrix), "x", len(mu_star_matrix[0]))

        return jsonify({
            'names': names,
            'output_vars': output_vars,
            'mu_star': mu_star_matrix,  # 矩阵
            'sigma': sigma_matrix  # 矩阵
        })
    except Exception as e:
        print("3PG模型分析错误:", str(e))
        return jsonify({"error": str(e)}), 400


# 3PG模型参数优化路由
@calc.route('/optimize_3pg', methods=['POST'])
def optimize_3pg():
    try:
        start_time = time.time()
        data = request.get_json()
        selected_params = data.get('selected_params', [])  # [{param_idx, var_idx}, ...]
        current_values = data.get('current_values', [])
        actual_results = data.get('actual_results', {})

        if not selected_params:
            return jsonify({"error": "没有选定需要优化的参数"}), 400
        if not current_values:
            return jsonify({"error": "当前参数值为空，无法进行优化"}), 400
        if not actual_results:
            return jsonify({"error": "未提供任何实际结果"}), 400

        # 验证索引
        for sel in selected_params:
            if sel['param_idx'] >= len(current_values) or sel['var_idx'] >= len(get_output_vars()):
                return jsonify({"error": f"索引超出范围: param_idx={sel['param_idx']}, var_idx={sel['var_idx']}"}), 400

        output_vars = get_output_vars()
        target_vars = [var for var in output_vars if var in actual_results]
        if not target_vars:
            return jsonify({"error": "实际结果中没有有效的输出变量"}), 400

        # 按输出变量分组选中的参数
        var_param_map = {}
        for sel in selected_params:
            var_idx = sel['var_idx']
            param_idx = sel['param_idx']
            if var_idx not in var_param_map:
                var_param_map[var_idx] = []
            var_param_map[var_idx].append(param_idx)

        # 构建搜索空间：针对所有选中参数
        unique_param_idxs = list(set(sel['param_idx'] for sel in selected_params))
        spaces = [Real(0, 1) for _ in unique_param_idxs]

        def objective(x):
            new_params = dict(default_parameters)
            for i, param in enumerate(current_values):
                new_params[param["name"]] = param["value"]
            # 更新优化参数
            for i, param_idx in enumerate(unique_param_idxs):
                param_config = available_params[param_idx]
                lb, ub = param_config["lower_bound"], param_config["upper_bound"]
                new_params[param_config["name"]] = lb + x[i] * (ub - lb)  # 反归一化

            model = create_3pg_model(params=new_params)
            model_results = model.run3PG(output_vars=target_vars)

            total_error = 0
            for var_idx, param_idxs in var_param_map.items():
                var = output_vars[var_idx]
                if var in actual_results:
                    if var in model_results and len(model_results[var]) > 0:
                        predicted = model_results[var][-1]
                        actual = actual_results[var]
                        total_error += (predicted - actual) ** 2
                    else:
                        total_error += actual_results[var] ** 2
            return total_error / len(target_vars)

        # 执行优化
        res = gp_minimize(objective, spaces, n_calls=50, n_random_starts=10, noise=1e-6, random_state=0)

        # 构建优化结果矩阵
        optimal_params = [[None for _ in range(len(output_vars))] for _ in range(len(current_values))]
        for i, param_idx in enumerate(unique_param_idxs):
            param_config = available_params[param_idx]
            lb, ub = param_config["lower_bound"], param_config["upper_bound"]
            optimized_value = round(lb + res.x[i] * (ub - lb), 6)
            for sel in selected_params:
                if sel['param_idx'] == param_idx:
                    optimal_params[param_idx][sel['var_idx']] = optimized_value

        end_time = time.time()
        elapsed_time = round(end_time - start_time, 2)

        return jsonify({
            'optimal_params': optimal_params,
            'optimal_value': round(res.fun, 6),
            'output_vars': output_vars,
            'elapsed_time': elapsed_time
        })
    except Exception as e:
        print("3PG模型优化错误:", str(e))
        return jsonify({"error": str(e)}), 400


@calc.route('/calculate', methods=['POST'])
def calculate():
    try:
        data = request.get_json()
        print("接收到计算请求：", data)  # 详细记录收到的数据
        parameters = data.get('parameters', [])

        # 使用用户调整的参数运行3PG模型
        user_params = dict(default_parameters)
        print("默认参数：", user_params)  # 打印默认参数

        # 清晰地记录每个参数的更新
        for param in parameters:
            param_name = param["name"]
            try:
                # 确保参数值是数值型
                param_value = float(param["value"])
                if param_name in user_params:
                    print(f"更新参数 {param_name}: {user_params[param_name]} -> {param_value}")
                    user_params[param_name] = param_value
            except (ValueError, TypeError) as e:
                print(f"参数值转换错误 {param_name}: {param['value']} - {str(e)}")

        print("更新后的参数：", user_params)  # 打印更新后的参数

        # 获取输出变量列表
        output_vars = get_output_vars()
        print("输出变量：", output_vars)

        # 创建3PG模型实例并使用用户参数
        model = create_3pg_model(params=user_params)

        # 运行模型并获取结果
        print("运行3PG模型...")
        model_results = model.run3PG(output_vars=output_vars)
        print("模型运行完成")

        # 处理结果
        final_results = process_model_results(model_results, output_vars)
        print(f"计算结果数量: {len(final_results)}")

        # 返回更详细的响应
        return jsonify({
            "model_results": final_results,
            "status": "success",
            "timestamp": datetime.now().isoformat()
        })
    except Exception as e:
        import traceback
        print(f"计算错误: {str(e)}")
        print(traceback.format_exc())  # 打印完整堆栈跟踪
        return jsonify({"error": str(e), "trace": traceback.format_exc()}), 400

# 进行最终结果的计算，展示灵敏度分析控制页面
@calc.route('/sa',methods=['POST','GET'])
def sa():
    if request.method == 'GET':
        benchmark_id = request.args.get('benchmark_id')
        project_id = request.args.get('project_id')
        project = ProjectsModel.query.get(project_id)

        nowproject = ProjectsModel.query.get(project_id)
        dataInput = {
            "Benchmark": benchmark_id
        }
        nowproject.data_input = json.dumps(dataInput)
        db.session.commit()

        params = load_parameters_config()

        # try:
        #     # 获取输出变量列表
        #     output_vars = get_output_vars()
        #     # 创建3PG模型实例
        #     model = create_3pg_model(Benchmark_id=benchmark_id)
        #     # 运行模型并获取结果
        #     model_results = model.run3PG(output_vars=output_vars)
        #     # 处理结果，只保留最后一个时间点的数据
        #     readable_results = process_model_results(model_results, output_vars)
        #     #print("3PG模型计算完成，结果数量:", len(readable_results))
        #
        #     return render_template('optimization.html',
        #                            available_parameters=params,
        #                            model_results=readable_results)
        # except Exception as e:
        #     print(f"3PG模型计算失败: {str(e)}")
        #     return render_template('optimization.html',
        #                            available_parameters=params,
        #                            model_results={})
        # 获取输出变量列表
        output_vars = get_output_vars()
        # 创建3PG模型实例
        model = create_3pg_model(Benchmark_id=benchmark_id)
        # model = create_3pg_model()
        # 运行模型并获取结果
        model_results = model.run3PG(output_vars=output_vars)
        # 处理结果，只保留最后一个时间点的数据
        readable_results = process_model_results(model_results, output_vars)
        # print("3PG模型计算完成，结果数量:", len(readable_results))

        return render_template('optimization_new.html',
                               available_parameters=params,
                               model_results=readable_results,
                               project=project)
    else:
        pass
